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Abstract. Studies of pairing correlations in ultrasmall metallic grains have 
commonly been based on a simple reduced BCS-model describing the scat- 
tering of pairs of electrons between discrete energy levels that come in time- 
reversed pairs. This model has an exact solution, worked out by Richardson 
in the context of nuclear physics in the 1960s. Here we give a tutorial intro- 
duction to his solution, and use it to check the quality of various previous 
treatments of this model. 



1. Introduction 

Recent experiments by Ralph, Black and Tinkham, involving the observa- 
tion of a spectroscopic gap indicative of pairing correlations in ultrasmall 
Al grains [1], have inspired a number of theoretical [2]-[ll] studies of how 
superconducting pairing correlations in such grains are affected by reduc- 
ing the grains' size, or equivalently by increasing its mean level spacing 
d oc Vol -1 until it exceeds the bulk gap A. In the earliest of these, a 
grand-canonical (g.c.) BCS approach [2, 3, 4] was applied to a reduced 
BCS Hamiltonian for uniformly spaced, spin-degenerate levels; it suggested 
that pairing correlations, as measured by the condensation energy , van- 
ish abruptly once d exceeds a critical level spacing d c that depends on the 
parity (0 or 1) of the number of electrons on the grain, being smaller for 
odd grains (df ~ 0.89A) than even grains (d^ ~ 3.6A). A series of more 
sophisticated canonical approaches (summarized in Section 3 below) con- 
firmed the parity dependence of pairing correlations, but established [6]- [11] 
that the abrupt vanishing of pairing correlations at d c is an artifact of g.c. 
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treatments: pairing correlations do persist, in the form of so-called fluctua- 
tions, to arbitrarily large level spacings, and the crossover between the bulk 
superconducting (SC) regime (d <C A) and the fluctuation-dominated (FD) 
regime {d 3> A) is completely smooth [10]. Nevertheless, these two regimes 
are qualitatively very different [9, 10]: the condensation energy, e.g., is an 
extensive function of volume in the former and almost intensive in the lat- 
ter, and pairing correlations are quite strongly localized around the Fermi 
energy ep, or more spread out in energy, respectively. 

After the appearance of all these works, we became aware that the re- 
duced BCS Hamiltonian on which they are based actually has an exact 
solution. It was published by R. W. Richardson in the context of nuclear 
physics (where it is known as the "picket-fence model"), in a series of pa- 
pers between 1963 and 1977 [12]-[20] which until very recently seem to have 
completely escaped the attention of the condensed matter community. In 
this work, we (i) give a tutorial introduction (with no pretense of rigor) 
to his solution, and (ii) compare the results of various previously-used ap- 
proximations against the benchmark set by the exact solution, in order 
to gauge their reliability for related problems for which no exact solutions 
exist [21, 22]. 

2. Richardson's Exact Solution 

2.1. REDUCED BCS MODEL 

Ultrasmall superconducting grains are commonly described [2]- [11] by a 
reduced BCS model, 

H = E A° c i° 4+ c l- c j- c j+' w 

jo ij 

for a set S of Ns pairs of time-reversed states \j, ±) labeled by a discrete 
index j = 1, . . . , Ns, with energies £j and coupling g = Xd, where d is the 
mean level spacing and A a dimensionless coupling constant. Unbeknownst 
to the authors that have studied this model recently, Richardson had long 
ago solved it exactly, for an arbitrary set of levels £j (degenerate levels 
are allowed, but are to be distinguished by distinct j-labels, i.e. they have 
Si = £j for i^j). 

The first step is to note that singly- occupied levels do not participate 
in the pairscattering described by H, and by the Pauli principle remain 
"blocked" [23] to such pairscattering; the labels of such levels are therefore 
good quantum numbers. A general eigenstate of H thus has the form 

\n,B) = Y[cl\y n )u, (2) 

ieB 
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U n 

\y n )u = E ^•••,i«)LK=ii°>- ( 3 ) 

ji,—,3n v=\ 

This describes N = In + b electrons, b of which sit in a set B of singly- 
occupied, blocked levels, thereby contributing 6b = J2ieB £ i t° the eigenen- 
ergy, while the remaining n pairs of electrons, created by the pair operators 
bj = c!j + Cj_, are distributed among the remaining set U = S\B of Njj = 

Ns — b unblocked levels, with wave function . . . , j n ) (Y^ = J2j^B 

denotes a sum over all unblocked levels). The dynamics of these pairs is 
governed by 

u 

^ = E( 2e A- 9)b\bj, (4) 

ij 

and writing the eigenenergy of \n, b) as 6 n + 6b, the state \^ n )u satisfies 

u 

H u \y n ) u = 6 n \V n ) u , J2 b ] b j\^n)u = n\V n )u . (5) 

3 

Diagonalizing Hjj would be trivial if the 6's were true bosons. However, they 
are not, and in the subspace spanned by the set U of all non-singly-occupied 
levels, instead satisfy the "hard-core boson" relations, 

bf = 0, [b h b),]=8 jj ,(l-2b%), [b%,b],} = 8 jr b], (6) 

which reflect the Pauli principle for the fermions they are constructed from. 
In particular, 6j = implies that only those terms in (3) are non-zero for 
which the indices ji, . . . j n are all distinct. 

In his original publications [12, 13, 14], Richardson derived a Schrodinger 
for . . . ,j n ) and showed that its exact solution was simply a general- 
ization of the form that ip(ji, ■ ■ ■ ,j n ) would have had if the 6's had been 
true (not hard-core) bosons. With the benefit of hindsight, we shall here 
follow an alternative, somewhat shorter root, also due to Richardson [24]: 
we first consider the related but much simpler case of true bosons and 
write down the generic form of its eigenstates; we then clarify why this 
form fails to produce eigenstates of the hard-core boson Hamiltonian; and 
having identified the reason for the failure, we show that (remarkably) only 
a slight generalization is needed to repair it and to obtain the sought-after 
hard-core-boson eigenstates. 

2.2. TRUE BOSONS 

Let bj denote a set of true bosons (i.e. = 5jj>), governed by a Hamil- 

tonian Hjj of precisely the form (4), with bj — > bj. This problem, being 
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quadratic, can be solved straightforwardly by any number of methods. The 
solution is as follows: Hjj can be written as 

Hu = Y, EjBjBj + const - ( 7 ) 
J 

where the new bosons Bj (with normalization constants Cj) are given by 
u l\ 1 u 1 

B ' J = 9Cj ^^Tj> WjY^^-Ej)^ (8) 

and the boson eigenenergies E j are the roots of the eigenvalue equation 

This is an equation of order Njj in Ej. It thus has Njj roots, so that the 
label J runs from 1 to Njj. As the coupling g is turned to 0, each Ej 
smoothly evolves to one of the bare eigenenergies £j. A general n-boson 
eigenstate of Hjj and its eigenenergy £ n thus have the form 

n n 

l*n)i/ = A b}Jo> , 4 = (io) 

where the n indices J±, . . . , J n that characterize this state need not all be 
distinct, since the B ] j are true bosons. 

2.3. COMPLICATIONS ARISING FOR HARD-CORE BOSONS 

Let us now return to the hard-core boson Hamiltonian Hjj. Its eigenstates 
will obviously not be identical to the true-boson eigenstates just discussed, 
since matters are changed considerably by the hard-core properties of bj. To 
find out exactly what changes they produce, it is very instructive to take 
an Ansatz for \^> n )jj similar to (10) (but suppressing the normalization 
constants and taking all J v to be distinct), namely 

n U 

\*n)u = n B i\°) > With B J=T, 2£ .1 E ■ (H) 

and to check explicitly whether or not it could be an eigenstate of Hjj, i.e. 
to check under what conditions (Hjj — £ n )\^ n )u would equal zero, where 
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£n = J2u Ej v ■ To this end, we commute Hjj to the right past all the B Jv 
operators in \^f n )u, using 




n n f / v — 1 \ In 

Hu,u s i =e n< Mj n b 

v=l J u =l ^ \n=l J \ij,=v+1 

To evaluate the commutators appearing here, we write Hjj as 

u u 
H U = Y / tejbfa - gB^Bo , where Bj = £ b) , 



(12) 



(13) 



and use the following relations: 



1 - 2&tfc 



j"3 



2sj - Ej ' 



(14) 



[Hu,B\\ = EjB\ + b\ 



" 1 ~ 2b}bj 

92 r isi-Ej 



(15) 



Inserting these into (12) and using Hu\0) = and £ n = J2u ^J v i we find 



H u \y n ) u = £ n \^ n ) u + Y^ 



n ( (v-\ 

+e n< 

v=\ I \ v =l 



E 



r 2 ^ - e j» 



2e, - £j„ 



4 n 4jio) 
v^=i(^) / 

n <)}io). as) 



Now, suppose we do the same calculation for true instead of hard-core 
bosons (i.e. run through the same steps, but place a ~ on Hjj, bj, Ej and 
£ n ). Then the second line of (16) would be absent (because the by>j terms 
in the second of Eqs. (6) and (14) and in (15) would be absent); and the 
first line of (16) would imply that (Hu — £ n )\^?n)u = provided that the 
term in square brackets vanishes, which is nothing but the condition that 
the Ej satisfy the the true-boson eigenvalue equation of (9)! In other words, 
we have just verified explicitly that all true-boson states of the form (10) 
are indeed eigenstates of Hjj, provided that the Ej satisfy (9). Moreover, 
we have identified the term in second line of (16) as the extra complication 
that arises for hard-core bosons. 
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2.4. THE CURE: A GENERALIZED EIGENVALUE EQUATION 
Fortunately, this extra complication is tractable: first, we note that 



U 2gBlb\b 



E 



2sj - E Jv 



u 

E 

3 



t 



2s -j - E Jv 2e j - Ej 



b\ -b 



Ej v - Ej 



(17) 



The rightmost expression follows via a partial fraction expansion, and re- 
markably, contains only £?t operators and no more b^bjs. This enables us 

to eliminate the b^bjs from the second line of (16), by rewriting it as fol- 
lows (we commute its term in square brackets to the right, using a relation 
similar to (12), but with the commutator (17) instead of [Hu,Bj ]): 




|0) 



(18) 



(The last line follows by renaming the dummy indices u <-> fj, in the second 
line.) Substituting (18) for the second line of (16), we conclude that {Hjj — 
£n)\^n)u will equal zero provided that 



u 

i-E 
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2e, - E Jv 



+ E 



2<? 



Ej 



Ej v 



0, 



for v = 1, 



,n 



(19) 



This consitutes a set of n coupled equations for the n parameters Ej x , . . . , 
Ej n , which may be thought of as self-consistently-determined pair ener- 
gies. Eq. (19) can be regarded as a generalization of the true-boson eigen- 
value equation (9), and was originally derived by Richardson by solving 
the Schrodinger equation for the wave-function . . . ,j n ) of (3). It is 

truly remarkable that the exact eigenstates of a complicated many-body 
problem can be constructed by such a simple generalization of the solution 
of a quadratic (i.e. non-interacting) true-boson Hamiltonian! 

Below we shall always assume the e/s to be all distinct. Then there ex- 
ists a simple relation between the bare pair energies 2ej and the solutions 
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of (19): as g is reduced to 0, it follows by inspection that each solution 
{Ej 1 , . . . , Ej n } reduces smoothly to a certain set of n bare pair energies, 
say {2ej 1 , . . . , 2£j n }- Correspondingly, the state \^ n )u = \Ji, ■ ■ ■ Jn)u of 
(11) reduces smoothly to the state . . . j n )u = II™=i ^jjO) ( U P t° a nor ~ 
malization factor not shown here). Thus there is a one-to-one correspon- 
dence between the set of all states {| Ji, . . . , J n )u} and the set of all states 
• • • jn)u}- Since the latter constitute a complete eigenbasis for the n- 
pair Hilbert space defined on the set of unblocked levels U, the former do 
too. 

2.5. GROUND STATE 

For a given set of blocked levels B, the lowest-lying of all states |n, B), say 
\n, B)q, is obtained by using that particular solution Ej x , . . . Ej n for which 
the total "pair energy" £ n takes its lowest possible value (as g is increased, 
some of the Ejs become complex; however, they always occur in complex 
conjugate pairs, so that £ n remains real [17]). 

The lowest-lying of all eigenstates with n pairs and b blocked levels, say 
\n,b)c with energy S^(n), is that \n, B)g for which the blocked levels in 
B are all as close as possible to ef, the Fermi energy of the uncorrelated 
iV-electron Fermi sea \Fn). The Ej v for the ground state \n,b)c coincide 
at g = with the lowest n energies 2ej (j = 1, . . . , n), and smoothly evolve 
toward lower values as g is turned on. This fact can be exploited during 
the numerical solution of (19), which can be simplified by first making 
some algebraic transformations, discussed in detail in [15], that render the 
equations less singular. 

2.6. GENERAL COMMENTS 

Since the exact solution provides us with wave functions, it is in principle 
straightforward to calculate arbitrary correlation functions. Some such cor- 
relators are discussed by Richardson in [16, 17], who showed that they can 
be expressed in terms of certain determinants that are most conveniently 
calculated numerically. Moreover, it is natural to ask whether in the bulk 
limit, the standard BCS results can be extracted from the exact solution. 
Indeed they can, as Richardson showed in [20], by interpreting the problem 
of solving (19) for the Ej v as a problem in two-dimensional electrostatics. 
Exploiting this analogy, he showed that in the bulk limit (N$ — > oo at 
at fixed Nsd), Eqs. (19) reduce to the well-known BCS gap equation and 
the BCS equation for the chemical potential, and the condensation energy 
£§{n) (defined in Eq. (20) below) to its BCS result, namely -A 2 /2d. 
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3. Comparison with Other Approaches 

We now apply the exact solution to check the quality of results previously 
obtained by various other methods. Most previous works [2, 3, 4, 6, 7, 8, 9, 
10] studied a half-filled band with fixed width 2ujd of uniformly-spaced lev- 
els (i.e. £j = j d), containing N = 2n + b electrons. Then the level spacing is 
d = 2uj£)/N and in the limit d — > the bulk gap is A = oj£> sinh(l/A) _1 . Fol- 
lowing [9], we take A = 0.224 throughout this paper. To study the SC/FD 
crossover, two types of quantities were typically calculated as functions 
of increasing d/A, which mimics decreasing grain size: the even and odd 
(6 = 0, 1) condensation energies 

EC{n)=£°{n)-{F N \H\F N )- (20) 

and a parity parameter introduced by Matveev and Larkin (ML) [6] to 
characterize the even-odd ground state energy difference, 

A ML H = S?(n) - [£§{n) + ££{n + l)]/2 . (21) 

Following the initial g.c. studies [2]- [6], the first canonical study was that of 
Mastellone, Falci and Fazio (MFF) [7], who used Lanczos exact diagonal- 
ization (with n < 12) and a scaling argument to probe the crossover regime. 
Berger and Halperin (BH) [8] showed that essentially the same results could 
be achieved with n < 6 by first reducing the bandwidth and renormaliz- 
ing A, thus significantly reducing the calculational effort involved. To access 
larger systems and fully recover the bulk limit, fixed-n projected variational 
BCS wavefunctions (PBCS) were used in [9] (for n < 600); significant im- 
provements over the latter results, in particular in the crossover regime, 
were subsequently achieved in [10] using the density matrix renormaliza- 
tion group (DMRG) (with n < 400). Finally, Dukelsky and Schuck [11] 
showed that a self-consistent RPA approach, that in principle can be ex- 
tended to finite temperatures, describes the f.d. regime rather well (though 
not as well as the DMRG). 

To check the quality of the above methods, we [21, 22] computed E^(n) 
and A ML (n) using Richardson's solution (Fig. 1). The exact results (a) 
quantitatively agree, for d — > 0, with the leading —A 2 /2d behavior for 
E^(n) obtained in the g.c. BCS approach [2, 3, 4], which in this sense is 
exact in the bulk limit, corrections being of order d°; (b) confirm that a 
completely smooth [10] crossover occurs around the scale d ~ A at which 
the g.c. BCS approach breaks down; (c) show that the PBCS crossover [9] 
is qualitatively correct, but not quantitatively, being somewhat too abrupt; 
(d) are reproduced remarkably well by the approaches of MFF [7] and BH 
[8]; (e) are fully reproduced by the DMRG of [10] with a relative error of 
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Figure 1. (a) The even and odd (6 = 0, 1) condensation energies of 
Eq. (20), calculated with BCS, PBCS and exact wave functions, as functions of 
d/A = 2sinh(l/A)/(2n + b), for A = 0.224. For comparison the dotted line gives the 
"bulk" result E^ ulk = -A 2 /(2d). (b) Comparison of the parity parameters A ML [6] of 
Eq. (21) obtained by various authors: ML's analytical result (dotted lines) [A(l — d/2A) 
for d <C A, and d/2 log(ad/A) for d ^> A, with a = 1.35 adjusted to give asymptotic 
agreement with the exact result]; grand-canonical BCS approach (dash-dotted line) [the 
naive perturbative result ^Xd is continued to the origin]; PBCS approach (short-dashed 
line); Richardson's exact solution (thick solid line); exact diagonalization and scaling by 
MFF (open circles) and BH (long-dashed line). 



< 10 4 for n < 400; our figures don't show DMRG curves, since they are 
indistinghuishable from the exact ones and are discussed in detail in [10]. 

4. Conclusions 

The main conclusion we can draw from these comparisons is that the 
two approaches based on renormalization group ideas work very well: the 
DMRG is essentially exact for this model, but the band-width rescaling 
method of BH also gives remarkably (though not quite as) good results 
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with rather less effort. In contrast, the PBCS approach is rather unreliable 
in the crossover region. 
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